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Abstract 

We  study  the  Rayleigh-Benard  stability  problem  for  a  fluid  confined  within  a  square  enclosure  subject  to  random  perturbations  in 
the  temperature  distribution  at  both  the  horizontal  walls.  These  temperature  perturbations  are  assumed  to  be  non-uniform  Gaussian 
random  processes  satisfying  a  prescribed  correlation  function.  By  using  the  Monte  Carlo  method  we  obtain  stochastic  bifurcation 
diagrams  for  the  Nusselt  number  near  the  classical  onset  of  convective  instability.  These  diagrams  show  that  random  perturbations 
render  the  bifurcation  process  to  convection  imperfect,  in  agreement  with  known  results.  In  particular,  the  pure  conduction  state 
does  no  longer  exist,  being  replaced  by  a  quasi-conduction  regime.  We  have  observed  subcritical  and  nearly  supercritical  quasi¬ 
conduction  stable  states  within  the  range  of  Rayleigh  numbers  Ra  =  0  -  4000.  This  suggests  that  random  perturbations  in  the 
temperature  distribution  at  the  horizontal  walls  of  the  cavity  can  extend  the  range  of  stability  of  quasi-conduction  states  beyond  the 
classical  bifurcation  point  Rac  =  2585.02.  Analysis  of  the  stochastic  bifurcation  diagrams  shows  the  presence  of  a  stochastic  drift 
phenomenon  in  the  heat  transfer  coefficient,  especially  in  the  transcritical  region.  Such  stochastic  drift  is  investigated  further  by 
means  of  a  sensitivity  analysis  based  on  functional  ANOVA  decomposition. 


1.  Introduction 

The  classical  stability  theory  of  Rayleigh-Benard  convec¬ 
tion  in  an  infinite  layer  of  fluid  confined  between  two  hori¬ 
zontal  isothermal  walls  with  constant  but  unequal  temperatures 
predicts  that  the  amplitude  of  the  motion  undergoes  a  bifurca¬ 
tion  as  the  Rayleigh  number  passes  through  the  critical  value 
Rac  =  1707.8  (see,  e.g.,  [1,  2]).  Such  bifurcation  characterizes 
the  transition  between  a  pure  conduction  state  and  convection. 
If  the  flow  is  laterally  confined  by  rigid  and  perfectly  insulating 
sidewalls  then  the  critical  Rayleigh  number  usually  increases 
[3,  4,  5,  6]  due  to  the  stabilizing  effects  of  the  finite  geometry. 
Furthermore,  if  there  is  small  heat  transfer  through  these  side- 
walls  so  that  the  boundary  conditions  are  inconsistent  with  a 
state  of  no-motion,  then  the  bifurcation  leading  to  convection 
is,  in  general,  replaced  by  a  smooth  transition  to  finite  ampli¬ 
tude  flow  [7].  Such  a  smooth  transition  has  been  also  predicted 
theoretically  for  thermal  convection  in  an  infinite  fluid  layer 
between  two  rigid  walls  with  different  mean  temperatures  and 
small  spatially  periodic  perturbations  [8],  Since  then,  a  consid¬ 
erable  research  effort  has  focused  on  examining  the  stability  of 
different  types  of  natural  convective  flows  subject  to  determin¬ 
istic  boundary  conditions  [9,  5,  10,  11,  12],  However,  not  as 
much  work  has  been  done  for  the  case  when  the  boundary  con¬ 
ditions  are  random  processes  of  finite  amplitude,  although  these 
results  would  bear  upon  the  importance  of  ignoring  uncertainty 
when  applying  classical  stability  results  in  real  situations,  both 
in  laboratory  experiments  and  elsewhere. 

Thus,  the  purpose  of  the  present  paper  is  to  examine  the 
effects  of  temperature  perturbations  on  the  classical  Rayleigh- 
Benard  stability  problem,  namely  an  unstably  stratified  fluid 
contained  between  two  smooth  horizontal  walls  with  different 


mean  temperatures.  In  particular,  we  will  study  the  prototype 
problem  of  a  square  enclosure  having  perfectly  insulating  lat¬ 
eral  sidewalls  and  determine  how  the  random  perturbations  in 
the  temperature  distributions  at  the  horizontal  walls  affect  the 
stability  and  the  branch  points  obtained  from  classical  bifur¬ 
cation  analysis.  Clearly,  when  no  variations  occur  along  the 
boundaries  convection  is  possible  only  when  the  Rayleigh  num¬ 
ber  is  greater  than  the  classical  critical  value  Rac  =  2585.02 
[6,  5,  13],  However,  when  random  temperature  variations  do 
occur  at  the  horizontal  walls,  the  bifurcation  process  leading 
to  convection  becomes  imperfect  [14]  and  the  subcritical  pme 
conduction  state  no  longer  exists,  being  replaced  by  a  quasi¬ 
conduction  regime  [8].  This  type  of  flow  is  characterized  by  a 
finite  -  though  perhaps  small  -  velocity  field  and  it  can  be  ob¬ 
served  even  at  low  values  of  Rayleigh  numbers. 

Many  important  questions  can  be  addressed  in  the  context 
of  stochastic  thermal  convection  driven  by  random  boundary 
conditions.  For  instance:  how  do  the  random  temperature  per¬ 
turbations  affect  stability  and  branch  points  obtained  from  clas¬ 
sical  bifurcation  analysis?  Is  there  any  connection  between  the 
stochastic  properties  of  the  temperature  perturbations  -  such  as 
correlation  length  and  amplitude  -  and  flow  stability?  Is  there 
a  preferential  correlation  length  enhancing  the  fluid  motion  and 
the  heat  transfer?  Is  it  possible  to  obtain  realizations  of  sta¬ 
ble  supercritical  quasi-conduction  states?  In  this  paper  we  will 
provide  an  answer  to  all  these  questions  by  employing  a  Monte 
Carlo  numerical  approach  [15,  16]  and  ANOVA  decomposition 
[17,18,19,20], 

This  paper  is  organized  as  follows.  In  section  2  we  formu¬ 
late  the  governing  equations  of  the  system,  i.e.,  the  Oberbeck- 
Boussinesq  approximation  to  convection  via  the  vorticity  trans- 
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where  <//(x,  v;  at)  and  T ( x ,  v;  at)  denote  the  streamfunction  and 
the  temperature  fields  while  Ra  and  Pr  are  the  Rayleigh  and 
the  Prandtl  numbers,  respectively.  The  variable  at  appearing 
in  i//(x,  y;  at)  and  T (x,  y;  at)  identifies  a  possible  outcome  of  the 
streamfunction  and  the  temperature  for  a  specific  realization  of 
the  random  temperature  distributions  at  the  horizontal  walls. 
All  the  quantities  have  been  made  dimensionless  by  scaling 
lengths  with  the  side  length  of  the  cavity  L,  streamfunction  with 
the  kinematic  viscosity  v,  time  with  L2/v  and  temperature  with 
a  reference  temperature  difference  A T>,  which  is  defined  to  be 
the  difference  between  the  spatial  averages  of  the  two  temper¬ 
ature  processes  at  the  horizontal  walls.  With  this  rescaling,  the 
Rayleigh  and  the  Prandtl  numbers  are  obtained  as 


Ra  = 


g/3L3ATr 


V 

Pr  = 

a 


(3) 


Figure  1 :  Schematic  of  the  dimensionless  geometry  and  dimensionless  temper¬ 
ature  boundary  conditions.  The  random  perturbations  gi  and  g2  are  assumed 
to  be  zero  mean  Gaussian  processes.  The  velocity  boundary  conditions  are  of 
no-slip  type,  i.e.,  iJj  =  difz/dx  =  dif//dy  =  0  at  the  solid  walls. 


port  equation  [21,  22],  In  section  3  we  characterize  the  random 
temperature  perturbations  at  the  horizontal  walls  of  the  cavity 
in  terms  of  a  Karhunen-Loeve  expansion  satisfying  a  prescribed 
Gaussian  correlation  function.  In  section  4  we  investigate  the 
effects  of  these  perturbations  -  parametrized  in  terms  of  their 
correlation  length  and  amplitude  -  on  the  onset  of  convective 
instability  and  we  determine  useful  stochastic  bifurcation  dia¬ 
grams  for  the  Nusselt  number  near  the  onset  of  convective  insta¬ 
bility.  The  existence  of  supercritical  quasi-conduction  states  is 
discussed  in  section  5.  By  using  the  ANOVA  method  in  section 
6  we  study  the  sensitivity  of  the  integrated  Nusselt  number  with 
respect  to  variations  in  the  amplitude  of  different  harmonics  ap¬ 
pearing  in  the  random  temperature  distributions  at  the  horizon¬ 
tal  walls.  This  allows  us  to  identify  the  most  effective  spatial 
frequency  enhancing  the  heat  transfer  coefficient.  Finally,  the 
main  findings  and  their  implications  are  summarized  in  section 
7.  We  also  include  two  brief  appendices  dealing  with  the  inte¬ 
gral  representation  of  the  Oberbeck-Boussinesq  equations  and 
the  description  of  the  ANOVA  technique  for  sensitivity  analy¬ 
sis,  respectively. 


where  g,  [3  and  a  are  the  acceleration  of  gravity,  the  isobaric 
compressibility  coefficient  and  the  thermal  diffusivity  of  the 
fluid,  respectively.  We  notice,  that  this  type  of  non-dimension- 
alization  is  not  effective  when  the  average  temperature  is  the 
same  along  the  boundaries.  Indeed,  in  this  case  the  reference 
temperature  difference  ATr  becomes  0  but  we  still  could  have 
convection  due  to  temperature  variations  at  the  boundary.  In 
figure  1  we  show  a  sketch  of  the  geometry  and  the  boundary 
conditions  associated  with  the  system  (l)-(2).  As  easily  seen, 
the  natural  convection  problem  we  are  examining  is  a  classical 
one,  i.e.,  an  incompressible  fluid  within  a  square  cavity  heated 
from  below  and  cooled  from  above.  The  sidewalls  of  the  cavity 
are  assumed  to  be  adiabatic  while  the  horizontal  walls  are  sub¬ 
ject  to  random  temperature  fluctuations  whose  rigorous  mathe¬ 
matical  definition  will  be  given  in  the  subsequent  section.  The 
velocity  boundary  conditions  are  assumed  to  be  of  no-slip  type, 
i.e.  di/s/dx  =  di/s/dy  —  0  at  solid  walls.  At  this  point  it  is  conve¬ 
nient  to  transform  the  non-homogeneous  temperature  boundary 
conditions  into  homogeneous  ones.  This  is  easily  achieved  by 
defining  the  new  field 


def  , 


T*(x,y;  w)  =  T(x,y,  at)  +  (>•  -  1)  (gi  (x;  to)  +  1)  -  yg2  (x;  to)  , 

(4) 

where  gi  (x;  at)  and  g2  (x;  to)  are  random  processes  satisfying 
adiabatic  boundary  conditions  at  x  =  0  and  x  =  1,  i.e. 


dgi 

dx 


=  0,  for  i  =  1,2. 


A  =0,1 


(5) 


2.  Governing  equations 


Let  us  consider  the  two-dimensional  steady  state  dimen¬ 
sionless  form  of  the  Oberbeck-Boussinesq  approximation  via 
the  vorticity  transport  equation  in  streamfunction-only  formu¬ 
lation 


d>ft  d(V2i/r)  dif/  d  (VV) 

dy  dx  dx  dy 

dif/  dT  _  dif/  dT 
dy  dx  dx  dy 


-PrV4if/  +  RaPr —  , 
dx 

V27\ 


(1) 

(2) 


Equation  (4)  can  be  inverted  to  give 

T  =  T*  +a-y)(gl  +  l)  +  yg2. 


From  Eq.  (6)  we  obtain 


dT 

dx 

dT 

~dy 

w2t 


dT*  |  (dg2_  dg\  \  +  dgt_ 
dx  '  \  dx  dx  J  dx 


+  (g2~g l)  -  1 , 
dy 


V2r +y 


d2g2 

dx2 


&gA  ^£1 

<9x2  /  <9x2 


(6) 

(7) 

(8) 
(9) 
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Finally,  a  substitution  of  Eqs.  (7),  (8)  and  (9)  into  Eqs.  (1)  and 
(2),  respectively,  yields  the  system 


dp  d(vV) 

dy  dx 


dp  3  (VV) 
dx  dy 


=  -PrW4p 


+RaPr 


(dgi 

_  dg  i 

\  dx 

dx 

dx  / 

(10) 


dip  /  dT*  |  (dg2 
dy  \  dx  +  '\  dx 


dgi\ 
dx  ) 


dgi  \  _  dp 
dx )  dx 

+V2r  +yl 


(£*<*-«■  >-) 
d2gi  _  d2g\  \  d2gi 
dx2  dx2  )  +  dx2 

(11) 


The  boundary  conditions  associated  with  Eqs.  (10)  and  (11)  are 
now  homogeneous.  In  Appendix  A  we  obtain  the  integral  rep¬ 
resentation  of  this  system  in  terms  of  eigenfunctions  of  proper 
eigenvalue  problems. 


3.  Characterization  of  temperature  perturbations  at  the 
horizontal  walls 


We  shall  assume  that  the  temperature  perturbations  gi  (x;  co) 
and  g 2  (x;  co)  are  zero  mean  random  processes  satisfying  adia¬ 
batic  boundary  conditions  at  x  =  0  and  x  —  1 .  In  order  to  rep¬ 
resent  these  processes  let  us  first  consider  a  suitable  orthonor¬ 
mal  basis  obtained  from  the  Sturm-Liouville  eigenvalue  prob¬ 
lem  [23,  24] 


d2p 

dx2 


+  a2p  =  0 , 


with 


dp  (0) 
dx 


dp(  1) 

dx 


(12) 


The  normalized  eigenfunctions  solving  (12)  are 


po(x)  =  1 ,  pn(x)  —  V2cos(«7rx)  n=  1,2,3,...  (13) 


Thus,  if  h  (x;  co)  is  a  zero-mean  process  in  [0, 1]  satisfying  adia¬ 
batic  boundary  conditions  at  x  =  0  and  x  =  1,  then  we  have  the 
following  spectral  representation 1  [25] 


h  (x;  co)  —  cr  z  ak  (oj)  (f)k  (x)  ,  (14) 

k=i 


where  cr  is  a  real  parameter  that  characterizes  the  amplitude  of 
the  process  while 


-if 

<r  Jo 


a/i  (co)  =  —  h  (x;  co)  pk  (x)  dx 


(15) 


are  uncorrelated  random  variables.  The  autocorrelation  of  the 
process  h  (x;  co)  has  the  obvious  representation 


def  </z(x;  co)h(x'\  co)) 
C  (x,  x  )  =  — 


crA 


'YJ(a2l)pn(x)pn(x')  ,  (16) 


n=  1 


Note  that  all  the  basis  functions  </>k(x)  (except  (k)  I  integrate  to  zero  and 
satisfy  adiabatic  boundary  conditions  at  x  =  0  and  x  =  1 . 


where  (•)  denotes  the  average  with  respect  to  the  joint  probabil¬ 
ity  measure  of  the  variables  \ak(co)}.  An  important  question  at 
this  point  is:  if  we  arbitrarily  prescribe  a  symmetric  autocorre¬ 
lation  function,  say  C*(x,y),  can  we  determine  a  set  of  uncor¬ 
related  random  variables  a*k  (co)  such  that  (16)  is  satisfied?  The 
answer  is  obviously  affirmative,  provided  the  prescribed  auto¬ 
correlation  satisfies  the  boundary  conditions 


dC*(x,y) 


dx 


=  0 , 


A  =0,1 


Vy  e  [0, 1] 


as  well  as  the  zero-mean  constraint 

-l 


Jo 


C*(x,y)dx  =  0  Vy  e  [0, 1] . 


(17) 


(18) 


If  C*(x,y)  does  not  satisfy  such  conditions  then  it  is  possible 
to  enforce  them  through  projection.  To  this  end,  let  us  first 
consider  the  (positive)  Fourier  coefficients 


<K) 


r*  1  r*  1 
def  f 

Jo  Jo 


C*  (x,  x')  pn  (x)  pn  (x')  dxdx  , 


n  >  1  (19) 


obtained  by  projecting  the  arbitrarily  prescribed  kernel  C*(x,  x') 
onto  the  basis  {p^}.  This  operation  basically  removes  every  spa¬ 
tial  gradient  at  the  boundaries  x  =  0  and  x  =  1  and  makes  the  as¬ 
signed  correlation  zero  spatial  mean,  in  the  sense  of  (18).  Next, 
let  us  consider  the  spectral  expansion  of  the  kernel  C*  (x,  x')  in 
terms  of  its  (positive)  eigenvalues  A ^  and  eigenfunctions  pk 


C*  (x,  x')  =  ^  Aki//k  (x)  \/fk  (x')  .  (20) 

k=  1 


A  substitution  of  this  expression  into  (19)  immediately  yields 


°°  r 

tt 


pk  (x)  pn  (x)  dx 


n  >  1  . 


(21) 


At  this  point  it  is  easy  to  check  that  if  {^-(w)}  is  any  set  of  zero- 
mean  and  uncorrelated  random  variables  with  unit  variance  (i.e. 
(£?)  =  1)  then  the  process 


h  (x;  a))  =  cr  ^<^>1/2&  (oj)  (pk  (x)  (22) 

k=  1 

satisfies  the  boundary  conditions  at  x  =  0  and  x  =  1  as  well  as 
the  zero  spatial  mean  constraint  and  it  has  the  following  corre¬ 
lation  function 


C*  (x,  x')  =  'Yj(bl)pn  (x)  pn  (x')  .  (23) 

n—l 


The  technique  just  discussed  can  be  considered  as  particular 
case  of  the  spectral  transformation  method  [26,  27]  where  an 
assigned  correlation  kernel  is  generated  by  assigning  the  spec¬ 
trum  relatively  to  a  specified  orthogonal  basis.  In  this  paper  we 
will  employ  the  following  Gaussian  correlation  function  (see 
[28]) 


C*  (x,  x')  =  exp 


(x-x')2 


(24) 
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(a)  (b) 


Figure  2:  (a)  Relative  energy  of  a  truncated  Karhunen-Loeve  expansion  of  the 
temperature  processes  at  the  horizontal  walls  as  a  function  of  the  number  of 
modes  retained  in  the  representation  for  different  correlation  lengths.  We  also 
show  the  relative  energy  cutoff  set  at  95%  (dashed  line),  (b)  Samples  of  temper¬ 
ature  perturbations  at  lower  and  upper  horizontal  walls  for  correlation  lengths 
lc  =  1  ( — ),  lc  =  0.5  (— )  and  lc  =  0.1  (-).  The  perturbation  amplitude  here  is 
set  at  5%  of  the  reference  temperature  difference. 


lc  oo  2  1 

0.5 

0.25 

0.1 

0.05 

0.025 

0.01 

M  12  3 

5 

9 

22 

44 

87 

199 

Table  1 :  Effects  of  correlation  length  on  the  dimensionality  of  the  temperature 
representation  at  each  horizontal  boundary.  The  energy  cutoff  is  set  at  95%  of 
the  total  energy  of  the  process. 

This  allows  us  to  represent  the  temperature  perturbations  g i  and 
g2  shown  in  figure  1  as 

M 

gi  (x;  oj)  =  a  (a>)  (pk(x)  ,  i  =  1,2  (25) 

k=l 

where  gr  (i  —  1,2;  k  —  1,2,  are  zero-mean  uncorre¬ 

lated  random  variables  with  unit  variance.  In  this  paper  we 
will  assume  that  are  standard  Gaussian  variables.  Several 
samples  of  the  processes  (25)  are  shown  in  figure  2  (b)  for  dif¬ 
ferent  correlation  lengths  /,  and  perturbation  amplitude  cr  set  at 
5%  of  reference  temperature  difference  between  the  horizontal 
walls.  Physically,  this  means,  e.g.,  a  temperature  perturbation 
with  amplitude  1  K  for  temperature  differences  of  about  20  K. 

We  remark  that  the  truncation  process  in  the  series  expansion 
(25)  has  to  be  performed  with  some  care,  in  such  a  way  that  the 
energy  of  the  neglected  modes  is  negligible.  To  this  end,  we 
examine  the  relative  energy  of  the  temperature  perturbations 

def  Ef  ( M )  def  9 

ef(M)=-f—,  where  Ef{M)=Y{b2n)  (26) 
Ef(° o) 

and  choose  the  total  number  of  terms  M  in  such  a  way  that  ef  is 
greater  than  a  specified  cutoff  value.  In  figure  2  (a)  we  show  the 
plots  of  ef  ( M )  corresponding  to  different  dimensionless  corre¬ 
lation  lengths  while  in  table  1  we  report  on  the  dimensionality 
M  -  i.e.  the  total  number  of  terms  -  of  the  spectral  representa¬ 
tion  (25)  for  a  95%  cutoff  threshold.  We  notice  that  as  the  cor¬ 
relation  length  goes  to  zero  the  temperature  perturbation  at  the 
boundaries  approaches  an  independent  increment  process  [29]. 
Even  for  temperature  perturbations  having  correlation  length  lc 
about  0.1  (scaled  on  the  side  length  of  the  cavity)  the  number 


of  expansion  terms  for  the  effective  representation  of  the  pro¬ 
cess  becomes  relatively  large.  Specifically,  since  we  have  two 
random  boundaries,  the  case  lc  =  0.1  results  in  a  stochastic 
system  forced  by  44  (22  +  22)  random  variables  (see  table  1). 
The  numerical  simulation  of  these  high-dimensional  problems 
requires  appropriate  techniques  [20,  30,  31,  32,  33,  34,  35,  36]. 
In  this  paper  we  will  employ  a  Monte  Carlo  method  but  also 
polynomial  chaos  with  adaptive  ANOVA  can  be  used  [19], 

4.  Stochastic  bifurcations  and  stability  of  steady  state  con¬ 
vection 

In  previous  work  [6]  we  have  obtained  bifurcation  diagrams 
for  natural  convective  flows  within  square  cavities  subject  to 
uniform  temperature  boundary  conditions.  We  have  observed 
the  coexistence  of  multiple  stable  steady  states,  in  agreement 
with  other  results  [12, 11,  10,  9],  for  the  same  values  of  Rayleigh 
and  Prandtl  numbers,  the  final  asymptotic  state  depending  on 
the  initial  flow  condition.  For  random  boundary  conditions, 
multiple  stable  states  can  still  exist  but  the  mechanism  of  their 
formation  is  substantially  different.  Indeed,  as  we  can  see  from 
Eqs.  (10)  and  (11),  the  random  perturbations  gi  and  g2  break 
the  symmetry  of  the  system,  i.e.,  the  steady-state  convection 
pattern  in  general  do  not  satisfy  the  discrete  symmetry  group 
described  in  [6],  Therefore,  the  symmetry-induced  multiplicity 
of  supercritical  states  in  this  case  is  replaced  by  a  more  physical 
ensemble  of  flows  in  a  one-to-one  correspondence  with  specific 
boundary  and  initial  conditions2.  We  have  identified  many  dif¬ 
ferent  steady-state  stable  convection  patterns  and  corresponding 
temperature  fields.  These  include  subcritical  and  supercritical 
quasi-conduction  states  for  which  the  kinetic  energy  of  the  flow 
turns  out  to  be  very  small.  In  figure  3  we  show  typical  temper¬ 
ature  fields  and  flow  patterns  corresponding  to  specific  realiza¬ 
tions  of  the  temperature  boundary  conditions.  The  rich  variety 
of  flows  associated  with  random  boundary  conditions  should 
be  compared  with  classical  results  of  convection  for  uniform 
temperature  distributions  (see,  e.g.,  [6,  5])  where  only  one  sub- 
critical  solution  (pure  conduction)  and  one  supercritical  solu¬ 
tion  (one-roll  pattern)  can  develop  within  the  range  of  Rayleigh 
numbers  considered  in  this  paper,  i.e.  Ra  -  0  -  4000. 

4.1.  Bifurcation  diagrams  for  the  Nusselt  number 

As  is  well  known,  a  sudden  change  in  the  slope  the  Nus¬ 
selt  number  versus  the  Rayleigh  number  usually  identifies  a 
transition  between  different  flow  states.  In  the  particular  case 
of  uniform  temperature  boundary  conditions  the  first  one  of 
these  transitions  characterizes  the  onset  of  convective  instabil¬ 
ity  [12,  5]  and,  for  the  geometry  shown  in  figure  1,  it  can  be 
clearly  identified  at  Rac  —  2585.02.  However,  in  the  presence 
of  random  temperature  perturbations  along  the  horizontal  walls 
of  the  cavity,  the  precise  determination  of  the  critical  Rayleigh 

“We  remark  that  for  very  specific  realizations  of  the  temperature  processes 
at  the  horizontal  walls,  convection  can  still  satisfy  the  discrete  symmetry  group 
described  in  [6].  However,  from  a  statistical  viewpoint  the  probability  that  this 
happens  is  zero. 
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Figure  3:  Typical  temperature  fields  (top  row)  and  streamlines  of  the  velocity  field  superimposed  to  the  modulus  of  velocity  (lower  row)  for  specific  realizations 
of  the  temperature  distribution  at  the  horizontal  walls  of  the  cavity.  The  amplitude  of  the  temperature  perturbations  is  set  at  5%  of  the  reference  temperature 
difference  while  the  correlation  length  is  lc  =0.1  (a),  lc  =  0.5  (b)  lc  =  0.25  (c).  Show  are:  (a)  a  subcritical  quasi-conduction  state  at  Ra  =  1946,  (b)  a  supercritical 
quasi-conduction  state  Ra  =  2650  and  (c)  a  fully  developed  one-roll  convection  pattern  at  Ra  =  3500. 


number  can  be  rather  difficult.  In  fact,  as  pointed  out  by  Ahlers 
et  al.  in  [14],  such  perturbations  render  the  bifurcation  pro¬ 
cess  to  convection  imperfect  and,  strictly  speaking,  a  critical 
Rayleigh  number  does  not  exist  in  the  usual  sense  since  con¬ 
vection  occurs  for  all  values  of  Ra.  However,  as  the  Rayleigh 
number  approaches  the  classical  critical  value,  the  amplitude  of 
convection  increases  greatly,  and  therefore  it  still  makes  sense 
to  define  a  “critical”  regime  near  the  classical  bifurcation  point. 

In  figure  4  we  show  the  bifurcation  diagrams  for  the  inte¬ 
grated  Nusselt  number 


lc  —  1  and  lc  =  0.5  are  very  similar.  This  can  be  explained 
by  noting  that  the  temperature  perturbations  at  the  horizontal 
walls  of  the  cavity  are  quite  similar  to  each  other  in  these  cases 
(see  figure  2  (b)).  Among  many  possible  convection  patterns, 
our  numerical  results  show  that  it  is  possible  to  obtain  realiza¬ 
tions  of  nearly  supercritical  (stable)  quasi-conduction  states. 
In  other  words,  it  seems  that  random  perturbations  can  stabi¬ 
lize  the  quasi-conduction  state  beyond  the  classical  bifurcation 
point.  This  rather  surprising  result  will  be  discussed  further  in 
the  next  section. 


dx  (27) 

y=0 

versus  the  Rayleigh  number.  These  diagrams  are  obtained  by 
first  sampling  the  temperature  distribution  at  the  horizontal  walls 
for  different  perturbation  amplitudes  and  correlation  lengths  and 
then  compute  the  corresponding  stable 3  convective  flow  through 
the  Galerkin  method  outlined  in  Appendix  A.  In  the  plots  of 
figure  4  we  also  include  the  classical  bifurcation  diagram  for 
deterministic  uniform  boundary  conditions  (dashed  lines).  This 
case  corresponds  to  lc  =  oo.  Note  that  the  bifurcation  diagrams 
obtained  for  temperature  perturbations  with  correlation  lengths 


...  .  def 

Nu  (to) 


.1  n 


dT(x,y;  to) 


dy 


3We  report  only  on  stable  steady  states.  Other  unstable  states  are  present  as 
well  but  these  are  not  shown  in  figure  4. 


4.2.  Statistical  analysis  of  the  heat  transfer 

In  figure  5  (a)  we  plot  the  probability  density  function  of  the 
integrated  Nusselt  number  at  Rayleigh  number  3000  (Pr  =  0.7) 
for  boundary  perturbations  with  different  correlation  lengths. 
Each  probability  density  is  estimated  through  a  non-parametric 
kernel  regression  method  based  on  the  available  temperature 
samples.  Specifically,  we  have  computed  105  flow  samples  at 
many  different  Rayleigh  numbers,  correlation  lengths  and  per¬ 
turbation  amplitudes  of  boundary  processes.  As  seen  from  fig¬ 
ure  5  (a),  random  temperature  perturbations  can  increase  or  de¬ 
crease  the  averaged  heat  transfer  relatively  to  the  uniform  case. 
In  a  mean  sense,  however,  it  turns  out  that  the  heat  transfer  is  en¬ 
hanced,  especially  in  the  transcritical  region  (see  figure  6).  Sim¬ 
ilarly,  in  figure  5  (b),  we  plot  the  probability  density  functions 
of  the  integrated  Nusselt  number  at  different  Rayleigh  numbers 
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Figure  4:  Bifurcation  diagrams  near  the  onset  of  convective  instability  for  random  temperature  perturbations  of  different  amplitudes  cr  and  correlation  lengths  lc.  Shown  are  the  integrated  Nusselt  numbers  versus  the 
Rayleigh  number  for  different  solution  samples.  The  dashed  line  in  each  plot  represents  the  classical  bifurcation  diagram  obtained  for  deterministic  uniform  temperature  conditions.  In  this  case  the  critical  Rayleigh 
number  is  Rac  =  2585.02. 


Figure  5:  (a)  Probability  density  functions  of  the  integrated  Nusselt  number  at 
Rayleigh  number  Ra  =  3000  for  boundary  perturbations  of  different  correlation 
lengths.  The  perturbation  amplitude  is  set  at  5%  of  the  reference  temperature 
difference.  The  vertical  line  indicates  the  deterministic  Nusselt  number  at  Ra  = 
3000  for  uniform  boundary  conditions,  (b)  Probability  density  functions  of  the 
integrated  Nusselt  number  at  different  Rayleigh  numbers:  Ra  =  1000  ( — ), 
Ra  =  2500  (•••),  Ra  =  3000  (-)  and  Ra  =  4000  (-  •  -).  The  temperature 
perturbations  have  correlation  length  lc  =  0.5  and  amplitude  cr  set  at  5%  of  the 
reference  temperature  difference. 

for  boundary  perturbations  with  correlation  length  lc  =  0.5  and 
amplitude  cr  set  at  5%  of  the  reference  temperature  difference. 
We  notice  that  at  Ra  -  1000  the  probability  density  of  Nu  is 
rather  peaked  around  Nu  =  1,  suggesting  a  high  probability 
of  quasi-conduction  regime.  In  the  transcritical  region  we  also 
observe  a  variation  of  the  probability  density  function  that  be¬ 
comes  approximately  Gaussian  when  convection  is  fully  devel¬ 
oped. 

Note  that  for  supercritical  flows,  the  probability  density  of 
the  integrated  Nusselt  number  is  continuously  supported.  This 
suggests  that  for  the  correlation  lengths  and  the  perturbation 
amplitudes  considered  in  this  paper  it  seems  that  there  exist 
only  one  possible  supercritical  convection  pattern,  i.e.  a  one- 
roll  flow.  In  other  words,  for  the  correlation  lengths,  the  per¬ 
turbation  amplitudes  and  the  range  of  Rayleigh  numbers  con¬ 
sidered  in  this  paper  the  ensemble  of  stable  flows  is  continuous 
and  composed  by  one -roll  patterns,  with  the  exception  of  some 
subcritical  quasi-conduction  states.  Next,  we  determine  the  av¬ 
erage  as  well  as  the  range  of  the  integrated  Nusselt  number  as  a 
function  of  the  Rayleigh  number.  This  study  helps  us  in  clarify¬ 
ing  if  the  correlation  length  of  the  temperature  perturbations  at 


Figure  6:  (a)  Mean  of  the  integrated  Nusselt  number  versus  the  Rayleigh  num¬ 
ber  for  boundary  perturbations  of  different  correlation  lengths.  The  perturba¬ 
tion  amplitude  is  set  at  5%  of  the  reference  temperature  difference  in  all  cases. 
We  see  that  random  temperature  perturbations  induce  a  stochastic  drift  phe¬ 
nomenon  in  the  transcritical  region,  (b)  Integrated  Nusselt  number  versus  the 
dimensionless  kinetic  energy  ( ec )  of  the  fluid  at  Prandtl  number  0.7.  The  am¬ 
plitude  of  the  boundary  perturbations  is  set  at  5%  and  the  correlation  length 
is  lc  =  1.  We  show  the  mean  (•••)  and  the  min-max  band  (-)  which  is 
parametrized  with  the  Rayleigh  number  Ra.  The  curves  at  constant  Ra  ( — ) 
are  simple  lines  due  to  the  very  high  correlation  coefficient  between  Nu  and  ec 
at  Prandtl  number  0.7. 

the  horizontal  walls  has  an  influence  on  the  averaged  heat  trans¬ 
fer  within  the  cavity.  To  this  end  we  examine  the  case  where  the 
perturbation  amplitude  is  set  at  5%  of  the  reference  temperature 
difference.  The  results  of  our  computations  are  shown  in  figure 
6  (a).  As  easily  seen,  random  temperature  perturbations  induce 
a  stochastic  drift  in  the  transcritical  region  yielding  to  an  in¬ 
crement  of  the  average  heat  flow.  This  increment  depends  on 
the  correlation  length  of  the  temperature  processes,  i.e.  there 
are  preferential  values  of  temperature  correlation  lengths  that 
trigger  convection  patterns  that  are  more  effective  for  what  con¬ 
cerns  the  heat  transfer.  Note,  however,  that  the  heat  transfer 
enhancement  is  rather  weak  in  all  cases  we  have  considered, 
quantifiable  in  approximately  10%  within  the  transcritical  re¬ 
gion.  Also,  when  convection  is  fully  developed  the  stochas¬ 
tic  drift  disappears  and  the  probability  density  of  the  integrated 
Nusselt  number  becomes  very  similar  to  a  Gaussian  distribution 
(see  figure  5  (b)). 

It  is  interesting  to  study  the  relation  between  the  integrated 
Nusselt  number  and  the  dimensionless  kinetic  energy  of  the 
fluid  in  more  detail.  Our  first  finding  is  that  the  correlation  coef- 
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ficient  between  these  two  quantities  is  approximately  one  in  all 
cases  we  have  considered  in  this  paper.  This  suggests  that  there 
exists  a  linear  relation  between  the  Nusselt  number  and  dimen¬ 
sionless  kinetic  energy  for  stochastic  convection  within  square 
cavities  at  fixed  Rayleigh  number  Prandtl  number  Pr  =  0.7. 
This  relation  is  shown  in  figure  6  (b)  where  we  plot  the  inte¬ 
grated  Nusselt  number  versus  the  kinetic  energy  of  the  fluid  for 
different  Rayleigh  numbers.  The  existence  of  a  linear  relation 
between  the  integrated  Nusselt  number  and  the  dimensionless 
kinetic  energy  implies  that  heat  transfer  is  primarily  determined 
by  advection,  even  in  the  quasi-conduction  regime. 

5.  Subcritical  and  supercritical  quasi-conduction  states 
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The  existence  of  subcritical  quasi-conduction  states  has  been 
theoretically  predicted  by  Kelly  &  Pal  in  [8]  for  an  infinite  layer 
of  fluid  with  small  periodic  temperature  variations  at  the  hor¬ 
izontal  walls.  By  means  of  perturbation  analysis,  they  have 
found  that  convection  can  occur  even  for  Rayleigh  numbers  less 
than  the  critical  one  (Rac  =  1707.8  for  the  infinite  layer).  The 
corresponding  Nusselt  number  in  this  case  is  a  function  of  the 
Rayleigh  number,  the  Prandtl  number  and  the  modulation  am¬ 
plitude.  The  perturbation  approach  of  Kelly  &  Pal,  however, 
cannot  be  easily  extended  to  the  present  flow  problem  because 
the  random  boundary  conditions  depend  on  many  variables  (see 
table  1)  and  it  is  not  easy  to  select  a  significant  perturbation  pa¬ 
rameter  quantifying  the  “amplitude  of  convection”4.  A  criterion 
to  identify  a  quasi-conduction  state  may  be  based  on  the  analy¬ 
sis  of  the  dimensionless  temperature  field  within  the  cavity.  In 
particular,  a  comparison  between  the  pure  conduction  solution 
and  the  convection  solution  can  reveal  if  there  is  a  significant 
temperature  transport  associated  with  the  fluid  motion.  The 
steady  state  pure  conduction  solution  can  be  easily  obtained  by 
integrating  the  Poisson’s  equation 


v2r  =  -y 


d2gi 

dx2 


d2g\ 
dx 2 


(28) 


with  homogeneous  boundary  conditions  (T*  =  0  at  the  horizon¬ 
tal  walls  and  dT* /dx  —  0  at  the  sidewalls  of  the  cavity).  Equa¬ 
tion  (28)  follows  from  Eq.  (2)  by  setting  i/j  =  0.  The  analytical 
solution  to  (28)  can  be  represented  in  terms  of  an  eigenfunction 
expansion  as 


T*  (x,y;  u> ) 


where 


4From  a  theoretical  viewpoint,  a  supercritical  stable  state  might  be  inves- 
tigated  by  analyzing  the  Oberbeck-Boussinesq  system,  in  any  representation. 
In  particular,  one  can  consider  the  integral  representation  obtained  in  appendix 
A,  expand  the  solution  near  =  0  and  bk  =  0  and  try  to  determine  whether 
there  exist  a  set  of  coefficients  for  which  the  real  part  of  the  largest  Jacobian 
eigenvalue  is  negative.  This  leads  to  a  complex  relation  between  the  forcing 
(buoyancy)  term  in  the  Navier-Stokes  equations  and  the  Rayleigh  number. 
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Figure  7:  Threshold  criterion  for  the  identification  of  quasi-conduction  states. 
These  diagrams  refer  to  the  case  where  the  boundary  perturbations  have  cor¬ 
relation  length  lc  =  0.5  and  perturbation  amplitude  set  at  5%  of  the  reference 
temperature  difference.  We  show  the  mean  Nusselt  number  (-),  the  minimum 
and  the  maximum  Nusselt  numbers  ( — )  and  the  classical  bifurcation  diagram 
(•  •  • )  obtained  for  spatially-uniform  deterministic  boundary  conditions.  Figure 
(b)  is  a  zoom-in  of  figure  (a). 


In  these  equations  y2  and  1 ),  denote  the  eigenvalues  and  the 
eigenfunctions  of  the  Helmholtz  operator,  respectively  (see  Ap¬ 
pendix  A.l  for  further  details).  Finally,  the  temperature  field 
associated  with  the  pure  conduction  state  corresponding  to  the 
perturbations  gi(x\  a>)  and  g2(x;  a>)  is  obtained  by  substituting 
Eq.  (29)  into  Eq.  (6).  A  simpler  criterion  to  identify  a  quasi¬ 
conduction  regime  is  based  on  the  integrated  Nusselt  number  it¬ 
self.  In  practice,  we  can  define  a  threshold  for  Nu  below  which 
we  can  state  that  convection  is  neglectable.  At  Prandtl  number 
0.7,  this  is  equivalent  to  selecting  a  threshold  for  the  dimen¬ 
sionless  kinetic  energy  of  the  fluid.  In  fact,  as  we  have  pointed 
out  in  the  previous  section,  the  integrated  Nusselt  number  and 
the  dimensionless  kinetic  energy  of  the  fluid  are  extremely  well 
correlated  in  all  cases  we  have  considered  in  this  paper.  The 
selection  of  a  threshold  value  for  the  integrated  Nusselt  num¬ 
ber,  however,  introduces  some  arbitrariness  in  the  definition  of 
quasi-conduction  states.  This  arbitrariness  is  of  the  same  type 
as  that  of  defining  a  critical  Rayleigh  number  in  the  presence  of 
random  boundary  conditions. 

Given  these  remarks,  let  us  set  the  threshold  Nu,r  —  1.02 


(a) 


Figure  8:  Probability  that  a  stable  quasi-conduction  state  develops  within  the 
cavity  as  a  function  of  the  Rayleigh  number  and  the  correlation  length  of  the 
temperature  processes  at  the  horizontal  walls  of  the  cavity.  The  perturbation 
amplitude  is  set  at  5%  of  the  reference  temperature  difference.  Figure  (b)  is  a 
zoom-in  of  figure  (a). 

for  quasi-conduction  states.  This  choice  is  based  on  a  careful 
analysis  of  the  temperature  fields  of  many  flow  samples  corre¬ 
sponding  to  different  realizations  of  boundary  conditions.  It  is 
clear  that  the  threshold  Nurr  -  1 .02  discriminates  among  those 
flows  with  heat  transfer  different  at  most  by  2%  with  respect 
to  pure  conduction.  In  figure  7  we  sketch  the  procedure  for 
the  identification  of  quasi  conduction-states  according  to  the 
proposed  criterion.  Clearly,  the  threshold  Nutr  =  1.02  is  also 
associated  with  a  quasi-conduction  kinetic  energy  band.  An 
analysis  of  the  stochastic  flow  field  near  the  onset  of  convec¬ 
tion  reveals  that  random  temperature  perturbations  at  the  hor¬ 
izontal  boundaries  can  stabilize  a  nearly  supercritical  quasi¬ 
conduction  regime.  This  region  is  indicated  in  figure  7  (b)  for 
boundary  perturbations  having  correlation  length  lc  =  0.5  and 
perturbation  amplitude  set  at  5%  of  the  reference  temperature 
difference.  Thus,  random  perturbations  basically  extend  the  do¬ 
main  of  stability  of  quasi-conduction  states  beyond  the  classi¬ 
cal  bifurcation  point.  The  flow  field  and  the  temperature  of  a 
supercritical  quasi-conduction  state  is  shown  in  figure  3  (b)  at 
Rayleigh  number  2650.  An  important  question  at  this  point 
is:  What  is  the  probability  that  a  supercritical  stable  quasi¬ 
conduction  state  develops  within  the  cavity? 

In  order  to  answer  this  question,  in  figure  8  we  plot  the 
probability  of  occurrence  of  quasi-conduction  states  within  the 
whole  range  of  Rayleigh  number  considered  in  this  paper,  for 
boundary  perturbations  of  different  correlation  lengths.  This 
probability  function  is  estimated  by  counting  the  relative  num¬ 


ber  of  quasi-conduction  states,  i.e.  the  states  whose  energy 
is  within  the  quasi-conduction  energy  band,  for  each  Rayleigh 
number.  As  easily  seen,  the  probability  curve  is  monotonic  and 
it  reaches  the  value  zero  (impossible  event)  approximately  at 
Ra  -  2800  in  all  cases.  We  also  notice  that  the  occurrence  of 
a  nearly  supercritical  quasi-conduction  state  is  rather  unlikely 
(see  figure  8  (b))  and  it  depends  on  the  correlation  length  of  the 
temperature  perturbations  at  the  horizontal  walls.  In  particu¬ 
lar,  smaller  correlation  lengths  yield  to  higher  probabilities  of 
supercritical  quasi-conduction  states.  Clearly,  all  these  results 
depend  on  the  choice  of  the  quasi-conduction  energy  band.  In 
other  words,  a  different  selection  of  the  threshold  for  the  Nusselt 
number  or  the  kinetic  energy  of  the  fluid  yields  quantitatively 
different  but  quantitatively  similar  conclusions. 

6.  Sensitivity  analysis 

In  this  section  we  employ  the  ANOVA  technique  [20,  17, 
37,  38,  32]  (see  also  appendix  B)  in  order  to  quantify  which 
harmonic  in  the  Fourier  series  representation  of  the  random 
temperature  boundary  conditions  enhances  the  heat  transfer  and 
triggers  the  transition  from  quasi-conduction  to  fully  developed 
convection  regimes.  This  sensitivity  study  allows  us  to  make 
inferences  about  the  most  important  unstable  modes  and,  in 
some  sense,  it  is  similar  to  the  perturbation  approach  adopted 
by  Kelly  &  Pal  [8]  for  the  infinite  fluid  layer. 

Let  us  consider  an  ANOVA  expansion  of  the  Nusselt  num¬ 
ber  in  terms  of  the  set  of  random  variables  representing  the  am¬ 
plitude  of  boundary  conditions 

2  M  2  M 

Nu  (O  =Nu0  +  Nui  (£)  +  Nuijdi,  (j) 
i=  1  i<j 

2  M 

+  2  &)  +  ■■■  ,  (31) 

i<j<k 

where 

(32) 

We  recall  that  M  depends  on  the  spatial  correlation  length  of 
the  temperature  process.  Specifically,  for  lc  -  0.5  -  which  is  the 
case  we  examine  here  -  we  obtain  a  total  number  of  10  (5  +  5) 
random  variables.  In  other  words,  the  nominal  dimension  [32] 
of  the  parameter  space  here  is  10  (see  table  1). 

The  sensitivity  (in  the  sense  of  Sobol  [39])  of  the  integrated 
Nusselt  number  with  respect  to  the  amplitude  of  the  boundary 
modes  can  be  studied  as  a  function  of  the  Rayleigh  number. 
This  provides  an  insight,  e.g.,  on  which  harmonic  of  the  temper¬ 
ature  distribution  at  the  boundaries  (first-order  interaction)  or 
combination  of  harmonics  (higher-order  interactions)  are  most 
important  in  the  transition  from  quasi-conduction  to  fully  de¬ 
veloped  convection.  The  results  of  this  study  are  summarized 
in  figure  9  where  we  plot  the  averaged  global  sensitivity  factors 
for  first-,  second-  and  third-order  interaction  terms  correspond¬ 
ing  to  all  five  parameters  [£j  \  ]  defining  the  random  tem¬ 

perature  process  at  the  lower  horizontal  wall.  These  sensitivity 
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factors  are  explicitly  defined  as 


(a) 


z(l)  def  CT2[NUj] 

1  ~  o~2[Nu]  ’ 

7(2)  def  (T2 

1  ~  cr2[Nu]  ’ 

^(3)  def  V  <j2tNuUk] 

'  “  Zj  o-2[Nu]  ’ 

J-fc 


(33) 

(34) 

(35) 


where  Nu„  Nujj  and  Nitijk  are  the  interaction  terms  in  Eq.  (31) 
while  cr2[-]  denotes  the  variance  operator  (see  Appendix  B  for 
further  details). 

As  seen  from  the  plots  of  Zyk)  ( k  =  1,2,3),  the  subcritical 
quasi-conduction  region  (Ra  <  Rac)  is  sensitive  to  variation  in 
the  amplitude  of  all  the  boundary  modes.  Beyond  the  classical 
bifurcation  point  we  also  see  that  there  is  transcritical  region 
(between  Ra  -  2700  and  Ra  -  3000)  which  is  very  sensitive 
to  variations  in  the  amplitude  of  the  first  boundary  mode,  i.e. 
cos(z*7rx).  We  recall  that  this  is  the  region  where  we  have  the 
stochastic  drift  phenomenon  for  the  Nusselt  number  (see  sec¬ 
tion  4.2).  Therefore,  we  can  conclude  that  the  average  heat 
transfer  enhancement  in  the  transcritical  region  is  due  to  a  low 
frequency  boundary  mode. 

Note  that  the  global  sensitivity  indices  associated  with  the 
first-order  interaction  terms  do  no  allow  us  to  identify  the  flow 
transition  in  a  clear  manner.  This  objective  is  achieved  indeed 
by  looking  at  the  higher-order  interaction  terms.  In  fact,  as  seen 
from  figure  9(b)  and  figure  9(c)  the  global  sensitivity  factors 
of  second-  and  third-order  interactions  terms  undergo  a  sud¬ 
den  jump  exactly  in  correspondence  with  the  classical  bifurca¬ 
tion  point.  This  suggests  that  the  interaction  between  different 
boundary  modes  is  switched  on  by  the  transition  and  the  result¬ 
ing  flow  becomes  rather  sensitive  to  variations  in  the  amplitude 
of  the  terms  associated  with  the  corresponding  harmonics.  We 
also  notice  that  there  is  a  bulk  phenomenon  in  the  sensitivity 
factors  within  the  region  of  fully  developed  convection,  i.e.  for 
Ra  >  3000.  This  suggests  that  in  such  region  the  Nusselt  num¬ 
ber  is  equally  sensitive  to  variations  in  the  amplitude  of  differ¬ 
ent  harmonics  of  the  temperature  expansion  at  the  lower  wall. 
This  is  expected  since  the  heat  transfer  in  the  fully  developed 
convection  region  in  primarily  determined  by  advection. 


7.  Summary 


Figure  9:  Averaged  global  sensitivity  indices  of  different  terms  in  the  ANOVA 
decomposition  of  the  Nusselt  number  for  variations  in  the  amplitude  of  the  har¬ 
monics  representing  the  temperature  boundary  condition  at  the  lower  horizontal 
wall.  Shown  are  sensitivities  of  (a)  first-order,  (b)  second-order  and  (c)  third- 
order  interaction  terms  versus  the  Rayleigh  number.  The  vertical  dashed  line 
identifies  the  classical  bifurcation  point  at  Rac  =  2585  .  It  is  seen  that  the  tran¬ 
sition  from  quasi-conduction  states  to  fully  developed  convection  is  captured 
by  the  second-  and  the  third-order  interaction  terms.  Also,  the  highest  sensitiv¬ 
ity  of  Nu  within  the  transcritical  region  is  achieved  by  the  variable  number  “1”. 
This  variable  characterizes  the  amplitude  of  the  lowest  frequency  mode  in  the 
Fourier  expansion  of  the  temperature  boundary  condition,  i.e.,  cos(7rx).  Thus, 
the  heat  transfer  enhancement  in  the  transcritical  region  (i.e.  the  stochastic  drift) 
is  mainly  influenced  by  such  harmonic. 


We  have  studied  the  Rayleigh-Benard  stability  problem  for 
fluid  confined  within  a  square  enclosure  subject  to  non-uniform 
random  perturbations  in  the  temperature  distribution  at  the  hor¬ 
izontal  walls.  These  temperature  perturbations  were  modelled 
as  Gaussian  processes  satisfying  a  Gaussian  correlation  func¬ 
tion.  We  have  simulated  the  Oberbeck-Boussinesq  equations 
and  computed  many  ensembles  of  realizations  of  the  natural 
convective  flow  within  the  cavity  by  sampling  the  temperature 
processes  at  the  boundaries  for  different  correlation  length  and 
amplitude.  This  allowed  us  to  obtain  stochastic  bifurcation  di¬ 
agrams  for  the  integrated  Nusselt  number  near  the  classical  on¬ 
set  of  convective  instability.  These  diagrams  show  that  random 
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perturbations  render  the  bifurcation  process  to  convection  im¬ 
perfect,  in  agreement  with  known  theoretical  results  [14],  In 
particular,  the  pure  conduction  state  does  no  longer  exist,  being 
replaced  by  a  quasi-conduction  regime.  We  have  observed  sub- 
critical  and  nearly  supercritical  quasi-conduction  stable  states 
within  the  range  of  Rayleigh  numbers  Ra  —  0  -  4000.  This 
suggests  that  random  temperature  perturbations  at  the  horizon¬ 
tal  walls  can  extend  the  range  of  stability  of  quasi-conduction 
states  beyond  the  classical  bifurcation  point.  However,  the  prob¬ 
ability  that  these  states  develop  within  the  cavity  is  rather  low. 
A  statistical  analysis  of  the  bifurcation  diagrams  near  the  clas¬ 
sical  onset  of  convection  shows  the  existence  of  a  stochastic 
drift  phenomenon  in  the  heat  transfer  coefficient,  especially  in 
the  transcritical  region.  The  increment  we  have  observed  in 
the  mean  Nusselt  number  is  about  10%  for  temperature  pertur¬ 
bations  having  a  correlation  length  comparable  with  the  side- 
length  of  the  cavity.  In  order  to  obtain  a  better  understanding  of 
this  phenomenon,  we  have  performed  a  sensitivity  analysis  of 
the  integrated  Nusselt  number  based  on  the  functional  ANOVA 
decomposition.  This  allowed  us  to  identify  which  harmonics  in 
the  random  temperature  distributions  at  the  horizontal  walls  are 
most  effective  in  enhancing  the  heat  transfer  coefficient.  The 
sensitivity  factors  corresponding  to  first-,  second-  and  third- 
order  interaction  terms  suggest  that  the  lowest-wavelength  har¬ 
monic  is  the  most  effective.  In  addition,  the  flow  transition  from 
quasi-conduction  to  fully  developed  convection  is  found  to  be 
accurately  captured  by  the  second-  and  the  third-order  interac¬ 
tion  terms. 
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A.  Integral  representations 

Let  us  consider  an  expansion  of  random  temperature  and 
velocity  fields  in  terms  of  normalized  eigenfunctions  t/r„  (x,  y) 
and  T,„  (x,y)5 

i//(x,y;a>)  =  ^an(aj)fin(x,y)  ,  (36) 

n=  1 
Nt 

T*  (x,y;  oj)  =  bm  (oj)  Tm  (x,  y)  .  (37) 

m=  1 

The  advantage  of  using  such  representations  is  that  they  au¬ 
tomatically  satisfy  all  the  boundary  conditions  as  well  as  the 
continuity  equation  [40,  21],  A  substitution  of  Eq.  (36)  and 
Eq.  (37)  into  ( 10)-(  11)  and  subsequent  Galerkin  projection  onto 
d/i;  and  l);,  respectively,  gives  the  following  system  of  alge¬ 
braic  equations  (repeated  indices  are  summed  unless  otherwise 
stated) 

ana,n'B„mk  PranC„k  RaPr(bnlDnk  +  A/)-)  0 ,  (38) 


5These  functions  are  obtained  in  Appendices  A.l  and  A. 2 


a I’nfinmk  @n  bj~nk  P nk )  Y /,  A  "b  A \k  0  ,  (39) 

where  the  coefficients  Nn ,  Air,  etc.,  are  defined  as 


.  ,  def 

Nk  = 


A  A  def 

Air  = 


Jf 


fd£2 

\  dx 

d2gi 
dx 2 


d-£f 


Jo  Jo 


\ipnij/kdxdy , 


Onk  = 


'0  JO 
■'  H  3?, 


r1  r 

Jo  Jo 


dx 


ipkdxdy , 


+ 


dgi 

dx 


if/kdxdy , 


Ykdxdy , 


V4t ljn\likdxdy , 
TP Ykdxdy , 

OX 


dgi'l  +  dM 

dx  )  dx  ) 


d^n 

dx 


(g2  -  g l) 


Ykdxdy , 


'  dij/„  3VVm 
^  dy  dx 

'  diKdY^  _ 
^  dy  dx 


dijjn  dW2t//„ 


dx  dy 
di An  dYm 


tf/kdxdy , 


dx  dy 


Ykdxdy . 


Also,  y?  denote  the  eigenvalues  of  the  Helmholtz  equation  (see 
Appendix  A.l).  The  nonlinear  system  (38)-(39)  can  be  solved 
numerically  with  accuracy  in  order  to  obtain  the  Fourier  coeffi¬ 
cients  a„(to)  and  bm((o)  for  each  realization  of  the  boundary  con¬ 
ditions.  Once  these  coefficients  are  available,  the  streamfunc- 
tion  and  the  temperature  fields  may  be  easily  recovered  from 
(36)  and  (37). 


A.l.  Temperature  expansion 

We  consider  an  eigenfunction  expansion  based  on  a  classi¬ 
cal  diffusion  problem  in  Cartesian  coordinates 

v2r  +  y2r  =  o  (40) 


with  homogeneous  boundary  conditions 


T*(x.  0)  =  T*(x,  1)  = 


3T*(0,y)  <9r*(l,y) 


=  0. 


(41) 


dx  dx 

A  separation  of  variables  in  Eq.  (40)  gives  the  following  two 
Sturm-Liouville  problems 


d2X  , 
d^+a~X~°' 


dX(  0)  dX(  1) 


dx 


dx 


=  0, 


d2Y 

— r-  +  b2Y  =  o ,  y (0)  =  y (i)  =  o , 

dy 


(42) 

(43) 


whose  solutions  are  the  well  known  [23,  24]  normalized  eigen¬ 
functions 


( 1  if  n  =  0 

Xn(x)  =  {  r 

|V2cos  {nnx)  if  n  =  1,2, 3,  ... 
Ym  (y)  =  V2  sin  (mny)  m  =  1,2,3,... 


(44) 

(45) 
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while  the  eigenvalues  are 


where  the  eigenvalues  a,-  are  solutions  of  the  transcendental 
equation 


a„  =  Tin  for  n  —  0, 1, 2, ... ,  (46) 

bm  =  Jim  for  m=  1,2,....  (47) 

This  implies  that  the  eigenvalues  of  (40)  are 

yL  =  n1  (»2  +  '»2)  •  (48) 

Thus,  the  two-dimensional  temperature  basis  function  can  be 
written  as 

T „  (x,  y)  =  Xi(n)  (x)  Y '](„)  (y) ,  (49) 

where  i(n)  and  j(n )  are  suitable  subsequences  obtained  accord¬ 
ing  an  ordering  of  y?.. 


A. 2.  Velocity  expansion 

All  the  velocity  boundary  conditions  are  of  no-slip  type. 
This  means  df/dx  =  dip/dy  =  0  everywhere  at  the  boundary. 
In  order  to  generate  a  divergence  free  basis  for  the  velocity  rep¬ 
resentation  satisfying  such  boundary  conditions  it  is  convenient 
to  consider  the  following  eigenvalue  problem  [40,  41,21,  12] 


d4if 

dx4 


dV 

dx4 


A  V, 


dili 

*Kx,  0)  =  /  =  o, 

ox 

dil/ 

my)  =  yr  =  0’ 

dy 


dil/ 

<Kx,  1)  =  f  -  =0, 
dx 

dil/ 

i«oj)  =  r  =  o- 

dy 


(50) 

(51) 

(52) 


This  eigenvalue  problem  is  symmetric6  and  separable.  A  sub¬ 
stitution  of  the  ansatz 


h  l a' )  -  I  _tan(Q''/2)  for  i  -  1, 3,5, ... 

tanh\T/  \  tan  (a,/2)  for  i  =  2,4,6,... 


(56) 


A  similar  solution  can  be  obtained  for  Y(y).  A  normalized  basis 
for  the  two-dimensional  streamfunction  can  be  obtained  as  a 
tensor  product  of  one-dimensional  bases  as 


<A n  (x,  y)  =  Xm  (x)  Y m  (y)  ,  (57) 


where  i(n)  and  j(n)  are  suitable  subsequences  obtained  accord¬ 
ing  to  an  eigenvalue  ordering 

<  =  4n)+0k»)-  (58) 


B.  ANOVA  decomposition  for  sensitivity  analysis 

The  key  idea  of  ANOVA  is  to  represent  a  high-dimensional 
function  f(x i ,  xi,  ■  ■  ■  ,  x,v)  in  terms  of  a  superimposition  of  func¬ 
tions  involving  a  lower  number  of  variables  (interaction  terms), 
and  then  truncate  the  series  at  specific  interaction  order.  Specif¬ 
ically,  the  ANOVA  expansion  of  an  A'-dimensional  scalar  func¬ 
tion  /  takes  the  form  [42] 


N  N 

fix  1,  x2,  ...,  XN)  =/o  +  ^  fi(Xi)  +  ^  fijiXj,  Xj ) 

*= 1  i<j 

N 

+  Yu  fjkixi,  Xj,  xk)  +  --  -  .  (59) 

i<j<k 


i//(x,y)  =  X(x)Y(y)  (53) 


into  (50)  yields  two  equivalent  eigenvalue  problems  for  X(x) 
and  Y{y)  in  the  form 


d4X 

dx4 


=  a4X, 


X(0)  =  X(l)  = 


dX(0)  dX(\) 


dx  dx 
The  general  integral  of  (54)  is  easily  found7  as 


=  0. 


(54) 

(55) 


X(x)  =  a  sin  (ax)  +  b  cos  (ax)  +  c  sinh  (ax)  +  d  cosh  (ax)  . 


The  function  /o  is  a  constant.  The  functions  fix,),  which  we 
shall  call  first-order  interactions,  give  the  overall  effects  of  the 
variables  x,  in  /  as  if  they  were  acting  independently  of  the 
other  input  variables.  The  functions  fijiXj,  xj)  describe  the  in¬ 
teraction  effects  of  the  variables  x,  and  x;-,  and  therefore  they 
will  be  called  second-order  interactions.  Similarly,  higher-order 
terms  reflect  the  cooperative  effects  of  an  increasing  number  of 
variables.  From  a  practical  viewpoint,  the  computation  of  the 
various  terms  in  the  ANOVA  expansion  can  be  performed  by  se¬ 
lecting  a  suitable  measure  space,  e.g.,  the  space  of  //-integrablc 
functions  in  the  hypercube  [0,  1  ],v,  where  p  denotes  an  integra¬ 
tion  measure.  In  this  case  we  have 


By  enforcing  the  boundary  conditions  (55) 
lowing  normalized  set  of  eigenfunctions 

we  obtain  the  fol- 

II 

O 

f1  fix  1... 
Jo 

.,  xt/)dp(x  i, ...,  xjv)  , 

(60) 

cos  [a,-  (x  -  1  /2)\  /  cos  [a,/2]  - 
cosh  [a,  (x  -  1/2)]  /  cosh  [a,/2] 

i  =  1,3,5,... 

fix,)  =  f  ••• 
Jo 

(V- 

Jo 

•5  Xr- 1  .  X/+J,  ...,  XN) 

Xj(x)  =  - 

sin  [a,-  (x  -  1/2)]  /  sin  [a, 72]  - 
sinh  [a,-  (x  -  1  /2)]  /  sinh  [a,/2] 

i  =  2,4,6,... 

dp(x i , . 

•  •  5  %i- 1  ?  %i+ 1 : 

,  ...,xjv)  -fo. 

(61) 

For  instance,  if  the  measure  p  is  selected  as 


6The  spectral  theory  for  linear  operators  in  a  Hilbert  space  guarantees  that 
the  eigenfunction  set  is  then  complete. 

7 It  is  sufficient  to  consider  a  solution  in  the  form  X(*)  =  eyx  to  obtain,  by 
substitution,  y4  =  a4,  i.e.,  y  =  {a,  —a,  ia ,  — ia } 


N 

dp(x i, ...,  x/y)  =  n  dxj  (62) 

i=l 


12 


then  we  obtain  the  classical  ANOVA-HDMR  method  [20].  Sim¬ 
ilarly,  if  we  set 


N 

dp(x i ,  •  •  •  ,  xN)  =  ]~ [  6(xj  -  c,)dxj ,  a  e  [0, 1]  (63) 

i=  1 

then  we  obtain  the  so-called  anchored  ANOVA  [43]  decompo¬ 
sition.  The  vector  (c i, ... ,cn )  in  this  case  is  known  as  anchor 
point  and  it  can  be  selected  according  to  many  different  crite¬ 
ria  (see,  e.g.,  [18]).  The  ANOVA  representation  of  a  field  can 
be  effectively  used  as  a  tool  for  sensitivity  analysis  [39,  44]. 
To  this  end,  let  us  first  recall  that  all  the  interaction  terms  (60), 
(61),  etc.,  in  the  ANOVA  expansion  (59)  are  mutually  orthogo¬ 
nal  with  respect  to  the  measure  p.  This  implies  that  the  variance 
of  /,  here  denote  as  cr2[/ ]  is  simply  the  sum  of  the  variances 
associated  with  each  interaction  term,  i.e.. 


N  N 

cr2[f]  -  X  ^Lfii  +  X  +  ■  • '  ’  (64) 

>'= 1  '<] 

where 

<r2[f]  =f  J  fdp  -  (J fdp J  .  (65) 

The  integrals  appearing  in  Eq.  (64)  can  be  computed  by  using  a 
multi-element  quadrature  formula  [32],  Following  Sobol  [39], 
we  shall  define  global  sensitivity  indices  as  the  ratio  between 
the  variance  of  each  term  in  the  ANOVA  decomposition  and 
the  total  variance  of  the  function  /,  i.e.. 


def  cr2  [/[]  def  O’  \fij] 

~  vHn  ’  u  - 


(66) 


From  Eq.  (64)  it  easily  follows  that 


N  N 

Xa’  -  XA’ (67> 

1=1  i<j 


Moreover,  we  shall  define  the  following  averaged  global  sensi¬ 
tivity  indices 


Z(1)  d=  R, 


(2)  def 


Z-J 

J 


IJ  ’ 


Z(3)  =f 


X  Ri>k  ’ 

j-k 


(68) 


representing  the  relative  importance  of  one  specific  parameter 
overall  the  others  at  a  prescribed  interaction  level.  With  the  aid 
of  Eqs.  (68)  we  can  study  which  input  variable  has  more  in¬ 
fluence  on  the  response  of  the  system.  For  instance,  we  can 
quantify  which  harmonic  in  the  Fourier  series  representation  of 
the  random  temperature  boundary  conditions  triggers  the  tran¬ 
sition  from  quasi-conduction  to  convection  or  affects  the  heat 
transfer  coefficient  to  the  greatest  extent. 
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